library("RColorBrewer")
library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
#library(SeuratWrappers)
library(harmony)
library(EnsDb.Hsapiens.v86)
library(stringr)
library(dplyr)
library(ggplot2)
library(patchwork)
library(kableExtra)
library(tidyverse)
set.seed(123)
# Paths
path_to_obj <- ("~/Documents/multiome_tonsil_Lucia/results/R_objects/14.tonsil_multiome_bcells_without_cluster4_doublets_normalized.rds")
path_to_markers<-("~/Documents/multiome_tonsil_Lucia/results/tables/tonsil_markers_bcell_withoutclust4_075.csv")
tonsil_wnn_bcell <- readRDS(path_to_obj)
tonsil_markers_075<-read_csv(path_to_markers)
## New names:
## * `` -> ...1
## Rows: 4815 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): ...1, gene
## dbl (6): p_val, avg_log2FC, pct.1, pct.2, p_val_adj, cluster
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
DimPlot(
tonsil_wnn_bcell,
group.by = "wsnn_res.0.075",
reduction = "wnn.umap",
pt.size = 0.1, label = T
)
Resolution 0.075
top5_tonsil_markers_075<-tonsil_markers_075 %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC)
top7_tonsil_markers_075<-tonsil_markers_075 %>% group_by(cluster) %>% top_n(n = 7, wt = avg_log2FC)
top10_tonsil_markers_075<-tonsil_markers_075 %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
df_top5<-as.data.frame(top5_tonsil_markers_075)
kbl(df_top5,caption = "Table of the top 5 marker of each cluster resolution 0.005") %>%
kable_paper("striped", full_width = F)
| …1 | p_val | avg_log2FC | pct.1 | pct.2 | p_val_adj | cluster | gene |
|---|---|---|---|---|---|---|---|
| COL19A1 | 0 | 2.3418794 | 0.799 | 0.221 | 0 | 0 | COL19A1 |
| STEAP1B | 0 | 2.1317977 | 0.521 | 0.126 | 0 | 0 | STEAP1B |
| CD69 | 0 | 1.6733680 | 0.519 | 0.174 | 0 | 0 | CD69 |
| CCSER1 | 0 | 1.5429622 | 0.646 | 0.303 | 0 | 0 | CCSER1 |
| LIX1-AS1 | 0 | 1.4711405 | 0.270 | 0.098 | 0 | 0 | LIX1-AS1 |
| ANK3 | 0 | 1.9694612 | 0.604 | 0.219 | 0 | 1 | ANK3 |
| SSPN | 0 | 1.7975224 | 0.446 | 0.099 | 0 | 1 | SSPN |
| AL355076.2 | 0 | 1.7115854 | 0.500 | 0.115 | 0 | 1 | AL355076.2 |
| HIPK2 | 0 | 1.6423040 | 0.379 | 0.144 | 0 | 1 | HIPK2 |
| PDE4D | 0 | 1.5809119 | 0.831 | 0.564 | 0 | 1 | PDE4D |
| HMGB2 | 0 | 3.0697315 | 0.956 | 0.135 | 0 | 2 | HMGB2 |
| TUBA1B | 0 | 2.9693211 | 0.969 | 0.222 | 0 | 2 | TUBA1B |
| H2AFZ | 0 | 2.7980292 | 0.968 | 0.233 | 0 | 2 | H2AFZ |
| TOP2A | 0 | 2.5647124 | 0.818 | 0.019 | 0 | 2 | TOP2A |
| HIST1H4C | 0 | 2.5492101 | 0.854 | 0.310 | 0 | 2 | HIST1H4C |
| MAML31 | 0 | 2.7180961 | 0.837 | 0.226 | 0 | 3 | MAML3 |
| AC023590.11 | 0 | 2.6018627 | 0.985 | 0.274 | 0 | 3 | AC023590.1 |
| AC104170.1 | 0 | 2.5051352 | 0.821 | 0.148 | 0 | 3 | AC104170.1 |
| RAPGEF51 | 0 | 2.3028431 | 0.933 | 0.233 | 0 | 3 | RAPGEF5 |
| LHFPL21 | 0 | 2.1961632 | 0.798 | 0.228 | 0 | 3 | LHFPL2 |
| IGHGP | 0 | 5.9845041 | 0.501 | 0.145 | 0 | 4 | IGHGP |
| IGHG1 | 0 | 5.8829304 | 0.638 | 0.291 | 0 | 4 | IGHG1 |
| IGHA1 | 0 | 6.4993828 | 0.646 | 0.424 | 0 | 4 | IGHA1 |
| IGKC | 0 | 5.9777857 | 0.961 | 0.915 | 0 | 4 | IGKC |
| IGLC2 | 0 | 5.9067244 | 0.894 | 0.755 | 0 | 4 | IGLC2 |
| MT-ND22 | 0 | 0.9424789 | 0.988 | 0.956 | 0 | 5 | MT-ND2 |
| ZBTB161 | 0 | 0.9873801 | 0.336 | 0.076 | 0 | 5 | ZBTB16 |
| AC105402.32 | 0 | 1.0221542 | 0.848 | 0.479 | 0 | 5 | AC105402.3 |
| GAB12 | 0 | 0.9870825 | 0.452 | 0.178 | 0 | 5 | GAB1 |
| FKBP51 | 0 | 0.9224334 | 0.572 | 0.335 | 0 | 5 | FKBP5 |
df_top7<-as.data.frame(top7_tonsil_markers_075)
df_mark<-as.data.frame(tonsil_markers_075)
kbl(df_top7,caption = "Table of the top 5 marker of each cluster resolution 0.005") %>%
kable_paper("striped", full_width = F)
| …1 | p_val | avg_log2FC | pct.1 | pct.2 | p_val_adj | cluster | gene |
|---|---|---|---|---|---|---|---|
| COL19A1 | 0 | 2.3418794 | 0.799 | 0.221 | 0 | 0 | COL19A1 |
| STEAP1B | 0 | 2.1317977 | 0.521 | 0.126 | 0 | 0 | STEAP1B |
| CD69 | 0 | 1.6733680 | 0.519 | 0.174 | 0 | 0 | CD69 |
| CCSER1 | 0 | 1.5429622 | 0.646 | 0.303 | 0 | 0 | CCSER1 |
| LIX1-AS1 | 0 | 1.4711405 | 0.270 | 0.098 | 0 | 0 | LIX1-AS1 |
| PTPRK | 0 | 1.4264085 | 0.471 | 0.130 | 0 | 0 | PTPRK |
| GAB1 | 0 | 1.3734829 | 0.305 | 0.108 | 0 | 0 | GAB1 |
| ANK3 | 0 | 1.9694612 | 0.604 | 0.219 | 0 | 1 | ANK3 |
| SSPN | 0 | 1.7975224 | 0.446 | 0.099 | 0 | 1 | SSPN |
| AL355076.2 | 0 | 1.7115854 | 0.500 | 0.115 | 0 | 1 | AL355076.2 |
| HIPK2 | 0 | 1.6423040 | 0.379 | 0.144 | 0 | 1 | HIPK2 |
| PDE4D | 0 | 1.5809119 | 0.831 | 0.564 | 0 | 1 | PDE4D |
| ZDHHC14 | 0 | 1.5635006 | 0.763 | 0.400 | 0 | 1 | ZDHHC14 |
| SYT1 | 0 | 1.5510644 | 0.248 | 0.034 | 0 | 1 | SYT1 |
| HMGB2 | 0 | 3.0697315 | 0.956 | 0.135 | 0 | 2 | HMGB2 |
| TUBA1B | 0 | 2.9693211 | 0.969 | 0.222 | 0 | 2 | TUBA1B |
| H2AFZ | 0 | 2.7980292 | 0.968 | 0.233 | 0 | 2 | H2AFZ |
| TOP2A | 0 | 2.5647124 | 0.818 | 0.019 | 0 | 2 | TOP2A |
| HIST1H4C | 0 | 2.5492101 | 0.854 | 0.310 | 0 | 2 | HIST1H4C |
| STMN1 | 0 | 2.5179082 | 0.941 | 0.090 | 0 | 2 | STMN1 |
| TUBB | 0 | 2.3594810 | 0.950 | 0.176 | 0 | 2 | TUBB |
| MAML31 | 0 | 2.7180961 | 0.837 | 0.226 | 0 | 3 | MAML3 |
| AC023590.11 | 0 | 2.6018627 | 0.985 | 0.274 | 0 | 3 | AC023590.1 |
| AC104170.1 | 0 | 2.5051352 | 0.821 | 0.148 | 0 | 3 | AC104170.1 |
| RAPGEF51 | 0 | 2.3028431 | 0.933 | 0.233 | 0 | 3 | RAPGEF5 |
| LHFPL21 | 0 | 2.1961632 | 0.798 | 0.228 | 0 | 3 | LHFPL2 |
| AFF21 | 0 | 2.0659813 | 0.927 | 0.255 | 0 | 3 | AFF2 |
| FGD61 | 0 | 2.0458983 | 0.860 | 0.226 | 0 | 3 | FGD6 |
| IGHGP | 0 | 5.9845041 | 0.501 | 0.145 | 0 | 4 | IGHGP |
| IGHG1 | 0 | 5.8829304 | 0.638 | 0.291 | 0 | 4 | IGHG1 |
| IGHG3 | 0 | 5.8542993 | 0.697 | 0.318 | 0 | 4 | IGHG3 |
| IGLC1 | 0 | 5.8772951 | 0.827 | 0.464 | 0 | 4 | IGLC1 |
| IGHA1 | 0 | 6.4993828 | 0.646 | 0.424 | 0 | 4 | IGHA1 |
| IGKC | 0 | 5.9777857 | 0.961 | 0.915 | 0 | 4 | IGKC |
| IGLC2 | 0 | 5.9067244 | 0.894 | 0.755 | 0 | 4 | IGLC2 |
| MT-ND22 | 0 | 0.9424789 | 0.988 | 0.956 | 0 | 5 | MT-ND2 |
| ZBTB161 | 0 | 0.9873801 | 0.336 | 0.076 | 0 | 5 | ZBTB16 |
| MT-ND12 | 0 | 0.9149171 | 0.988 | 0.966 | 0 | 5 | MT-ND1 |
| AC105402.32 | 0 | 1.0221542 | 0.848 | 0.479 | 0 | 5 | AC105402.3 |
| GAB12 | 0 | 0.9870825 | 0.452 | 0.178 | 0 | 5 | GAB1 |
| ST6GALNAC31 | 0 | 0.8516123 | 0.544 | 0.234 | 0 | 5 | ST6GALNAC3 |
| FKBP51 | 0 | 0.9224334 | 0.572 | 0.335 | 0 | 5 | FKBP5 |
#install.packages("htmlwidgets", type = "binary")
#install.packages("DT", type = "binary")
DT::datatable(df_top7)
DT::datatable(df_mark)
markers<-c("BANK1", "FCER2","CD3D", "IL7R","MKI67", "TOP2A","MARCKSL1", "RGS13", "LMO2", "CCDC88A","GNLY", "NKG7", "GZMK", "CD8A", "FCRL4", "FCRL5", "PLAC8", "SOX5","IGHG1", "IGHA1", "JCHAIN", "XBP1","LYZ", "S100A8","MT2A", "CD3D", "TRAC", "PCNA","CD19","CR2","MS4A1","RALGPS2","CD79A")
markerGenes <- unique(markers)
geneSym <- ifelse(test = !grepl('NA', markerGenes),
yes = sub('.*?-', '', markerGenes),
no = sub('-.*', '', markerGenes))
Idents(tonsil_wnn_bcell) <- "wsnn_res.0.075"
dot <- DotPlot(tonsil_wnn_bcell, features = unique(markers),cols = c("lightgrey", "blue"), cluster.idents = T) + theme(axis.text.x = element_text( size = 10, vjust = 0.8, hjust = 0.8)) + scale_x_discrete(labels= geneSym)+ggtitle("Known markers")
dot +
coord_flip() +
theme(axis.title = element_blank(), axis.text.y = element_text(size = 5))
Idents(tonsil_wnn_bcell) <- "wsnn_res.0.075"
dot.10 <- DotPlot(tonsil_wnn_bcell, features = unique(top10_tonsil_markers_075$gene),cols = c("lightgrey", "blue"), cluster.idents = T) + theme(axis.text.x = element_text( size = 10, vjust = 0.8, hjust = 0.8)) +ggtitle("res 0.1 top 10 of each cluster")
dot.5 <- DotPlot(tonsil_wnn_bcell, features = unique(top5_tonsil_markers_075$gene),cols = c("lightgrey", "blue"), cluster.idents = T) + theme(axis.text.x = element_text( size = 10, vjust = 0.8, hjust = 0.8))+ggtitle("res 0.1 top 5 of each cluster")
dot.10 +
coord_flip() +
theme(axis.title = element_blank(), axis.text.y = element_text(size = 5))
dot.5 +
coord_flip() +
theme(axis.title = element_blank(), axis.text.y = element_text(size = 7))
top7mark_cluster0<-top7_tonsil_markers_075[["gene"]][1:7]
top7mark_cluster1<-top7_tonsil_markers_075[["gene"]][8:14]
top7mark_cluster2<-top7_tonsil_markers_075[["gene"]][15:21]
top7mark_cluster3<-top7_tonsil_markers_075[["gene"]][22:28]
top7mark_cluster4<-top7_tonsil_markers_075[["gene"]][29:35]
top7mark_cluster5<-top7_tonsil_markers_075[["gene"]][36:42]
markers_gg <- function(x){purrr::map(x, function(x) {
p <- FeaturePlot(
tonsil_wnn_bcell,
features = x,
reduction = "wnn.umap",
pt.size = 0.1
)
p
})}
m<-c("PRDM1","XBP1","IRF4","MEF2B","BCL6")
DZ<-c("SUGCT", "CXCR4", "AICDA")
LZ<- c("CD83","BCL2A1")
GC<- c("MEF2B", "BCL6","IRF4")
PC<- c("PRDM1","SLAMF7", "MZB1", "FKBP11")
markers_gg(DZ)
## [[1]]
##
## [[2]]
##
## [[3]]
markers_gg(LZ)
## [[1]]
##
## [[2]]
markers_gg(GC)
## [[1]]
##
## [[2]]
##
## [[3]]
markers_gg(PC)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
## Pauli markers
naive_mem_bcell<-c("BANK1", "FCER2")
cd4_tcell<-c("CD3D", "IL7R")
dz_gc_bcell<-c("MKI67", "TOP2A")
lz_gc_bcell<-c("MARCKSL1", "RGS13", "LMO2", "CCDC88A")
cytotoxic<-c("GNLY", "NKG7", "GZMK", "CD8A")
memory_bcell<-c( "FCRL4", "FCRL5", "PLAC8", "SOX5")
pc<-c("IGHG1", "IGHA1", "JCHAIN", "XBP1")
myeloid<-c("LYZ", "S100A8")
poor_q_doublets <-c("FDCSP", "CLU", "CXCL13", "CR2")
doublet_proliferative_tcell<-c("MT2A", "CD3D", "TRAC", "PCNA")
Unk<-c("PTPRCAP", "CD37", "CD74")
PDC<-c("PTCRA", "LILRA4", "IRF7")
TCL1A:memory naive b cells IGHD:Group enriched (naive B-cell, memory B-cell) FCER2: (naive B-cell, memory B-cell) IGHM (naive B-cell, memory B-cell) IL4R:naive B-cell
markers_gg(c("TCL1A","IGHD","FCER2","IGHM","IL4R"))
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
markers_gg(naive_mem_bcell)
## [[1]]
##
## [[2]]
markers_gg(memory_bcell)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg(dz_gc_bcell)
## [[1]]
##
## [[2]]
markers_gg(lz_gc_bcell)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg(pc)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg(poor_q_doublets)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
Proliferating cell nuclear antigen (PCNA)
markers_gg(Unk)
## [[1]]
##
## [[2]]
##
## [[3]]
markers_gg(PDC)
## [[1]]
##
## [[2]]
##
## [[3]]
markers_gg( "MYO1E")
## [[1]]
markers_gg("COL4A3")
## [[1]]
FCRL4 is an immunoregulatory receptor that belongs to the Fc receptor-like (FCRL) family. In healthy individuals, FCRL4 is specifically expressed by memory B cells (MBCs) localized in sub-epithelial regions of lymphoid tissues. Expansion of FCRL4+ B cells has been observed in blood and other tissues in various infectious and autoimmune disorders ptprk: naive-memory
markers_gg(top7mark_cluster0)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
naive_markers<-c("CD79A", "CD79B", "BLNK")
memory_markers<-c("CD27")
markers_gg(naive_markers)
## [[1]]
##
## [[2]]
##
## [[3]]
markers_gg(memory_markers)
## [[1]]
markers_gg(c("MS4A1","NT5E"))
## [[1]]
##
## [[2]]
MS4A1:NAIVE-MEMORY B CELL
NT5E: NAIVE B CELL
markers_gg(top7mark_cluster1)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
markers_gg(top7mark_cluster2)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
markers_gg(top7mark_cluster3)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
INPP4B: Immune cell enhanced (memory CD4 T-cell)
AOAH: NK GNLY:NK
markers_gg(top7mark_cluster4)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
markers_gg(top7mark_cluster5)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
p <- FeaturePlot(
tonsil_wnn_bcell,
features = c("NFKB1","FOXO1"),
reduction = "wnn.umap",
pt.size = 0.1
)
p
FeaturePlot(
tonsil_wnn_bcell,
features = "percent.mt",
reduction = "wnn.umap",
pt.size = 0.1
)
VlnPlot(tonsil_wnn_bcell, features = "percent.mt", group.by = "wsnn_res.0.075", pt.size=0)
VlnPlot(tonsil_wnn_bcell, features = "percent_ribo", group.by = "wsnn_res.0.075", pt.size=0)
weight.vp<- VlnPlot(tonsil_wnn_bcell, features = "RNA.weight", group.by = "wsnn_res.0.075",pt.size = 0)
p.weight.vp<- VlnPlot(tonsil_wnn_bcell, features = "peaks.weight", group.by = "wsnn_res.0.075",pt.size = 0)
weight.vp + ggtitle("RNA modality weight")
p.weight.vp+ ggtitle("ATAC modality weight")
Idents(tonsil_wnn_bcell) <- "wsnn_res.0.075"
cell.num <- table(Idents(tonsil_wnn_bcell))
cell.num
##
## 0 1 2 3 4 5
## 15538 11029 7224 6704 1865 250
new.cluster.ids <- c("MBC", "NBC","GC/DZ", "GC/LZ", "GC/PC", "PC")
names(new.cluster.ids) <- levels(tonsil_wnn_bcell)
tonsil_wnn_bcell <- RenameIdents(tonsil_wnn_bcell, new.cluster.ids)
DimPlot(tonsil_wnn_bcell, reduction = "wnn.umap", label = TRUE, pt.size = 0.5)
MARKERS
Immature B cells express CD19, CD 20, CD34, CD38, and CD45R, T-cell receptor/CD3 complex (TCR/CD3 complex)
canonical_bcell_markers <-c("CD34", "CD38", "CD19")
monocytes_markers<-c("LYZ","S100A8")
naive_markers<-c("CD79A", "CD79B", "BLNK")
bib_Bcell_markers<-c("CD19","CR2","MS4A1","RALGPS2","CD79A")
bib_Tcell_markers<-c("CD3E","CD4","CD8A","FOXP3","IL17A")
markers_gg(naive_markers)
## [[1]]
##
## [[2]]
##
## [[3]]
markers_gg(bib_Bcell_markers)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
CD8+ T cell markers:“CD3D”, “CD8A” NK cell markers:“GNLY”, “NKG7”
markers_gg(bib_Tcell_markers)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
DimPlot(tonsil_wnn_bcell,
reduction = "wnn.umap",
pt.size = 0.1, label = T, group.by = "age_group")
DimPlot(tonsil_wnn_bcell,
reduction = "wnn.umap",
pt.size = 0.1, label = T, group.by = "Phase")
DimPlot(tonsil_wnn_bcell,
reduction = "wnn.umap",
pt.size = 0.1, label = T)